#ifndef FUNCTION_HPP
#define FUNCTION_HPP

#include <iostream>
#include <math.h>
#include <vector>
#include "norm.hpp"
using namespace std;

template <typename T>
vector<T> operator+(const vector<T> &lhs, const vector<T> &rhs)
{
    vector<T> res(lhs.size());
    for (size_t i = 0; i < res.size(); i++)
    {
        res[i] = lhs[i] + rhs[i];
    }
    return res;
}

template <typename T>
vector<T> operator-(const vector<T> &lhs, const vector<T> &rhs)
{
    vector<T> res(lhs.size());
    for (size_t i = 0; i < res.size(); i++)
    {
        res[i] = lhs[i] - rhs[i];
    }
    return res;
}

template <typename T>
vector<T> operator*(const T &lhs, const vector<T> &rhs)
{
    vector<T> res(rhs.size());
    for (size_t i = 0; i < res.size(); i++)
    {
        res[i] = lhs * rhs[i];
    }
    return res;
}

template <typename T>
ostream &operator<<(ostream &os, const std::vector<T> &v)
{
    os << v[0] << " " << v[1] << ";";
    return os;
}

vector<double> f(const vector<double> &v, double t)
{
    double mu = 0.012277471;
    vector<double> vec;
    vec.push_back(v[3]);
    vec.push_back(v[4]);
    vec.push_back(v[5]);
    vec.push_back(2 * v[4] + v[0] - (mu * (v[0] + mu - 1)) / pow((pow(v[1], 2) + pow(v[2], 2) + pow(v[0] + mu - 1, 2)), 1.5) - ((1 - mu) * (v[0] + mu)) / pow((pow(v[1], 2) + pow(v[2], 2) + pow(v[0] + mu, 2)), 1.5));
    vec.push_back(-2 * v[3] + v[1] - (mu * v[1]) / pow((pow(v[1], 2) + pow(v[2], 2) + pow(v[0] + mu - 1, 2)), 1.5) - ((1 - mu) * v[1]) / pow((pow(v[1], 2) + pow(v[2], 2) + pow(v[0] + mu, 2)), 1.5));
    vec.push_back(-(mu * v[2]) / pow((pow(v[1], 2) + pow(v[2], 2) + pow(v[0] + mu - 1, 2)), 1.5) - ((1 - mu) * v[2]) / pow((pow(v[1], 2) + pow(v[2], 2) + pow(v[0] + mu, 2)), 1.5));
    return vec;
}


class TimeIntegrator
{
public:
    TimeIntegrator(){}
    virtual void solver(){}
    virtual void error(){}
};

#endif